In [ ]:
sigma_clip
The following reduction code works directly from the full set of images in a given nights observations. As observations are added during the evening, the ImageFileCollection object defined at the beginning needs to be refreshed.
Files can be rsync'd from the whtobs computer using the following command:
$ rsync -av whtobs@taurus.ing.iac.es:/obsdata/whta/20170624/ 20170624/
where 20170624/ is set to the relevant observation night.
The starting point for this code was the very helpful scripts of Steve Crawford. However, various changes have been made and key steps (e.g. wavelength calibration/flat processing) differ.
Propagate uncertainties to create variance map
Flux calibration
In [2]:
import numpy as np
from multiprocessing import Pool
from functools import partial
import astropy.units as u
from astropy.io import fits
from astropy.modeling import models, fitting
from astropy.stats import sigma_clip, mad_std
from scipy.ndimage import binary_dilation
from astropy.utils.console import ProgressBar
import ccdproc
from ccdproc import ImageFileCollection, CCDData
def fit_chebyshev(row, degree=5, grow=3):
"""
Fit Chebyshev1D model to a CCD row, including masking of outlier pixels
Params
------
row : array,
Input CCD row to fit polynomial to.
degree : int,
Chebyshev Polynomial Order
grow : int,
Number of iterations to dilate around masked pixels
"""
fitter = fitting.LinearLSQFitter()
input_mask = row.mask
clipped = sigma_clip(row, stdfunc=mad_std)
clipped_pixels = np.array(clipped.mask+row.mask).astype('float')
clipped_pixels = binary_dilation(clipped_pixels, iterations=grow)
row[clipped_pixels==1] = np.median(row)
masked_row = np.ma.array(data=row,
mask=(clipped_pixels == 1),
fill_value=np.median(row))
x = np.arange(len(row))
model = models.Chebyshev1D(degree=degree)
fitted_model = fitter(model, x, row)
return fitted_model(np.arange(len(row)))
def fit_background_old(data, degree=5, grow=3, verbose=True):
"""
Background estimation for longslit CCD image
Params
------
ccd : array,
Input CCD data for background estimation.
degree : int,
Chebyshev Polynomial Order
grow : int,
Number of iterations to dilate around masked pixels
verbose : bool, default=True
If true, print progress bar
"""
fitted_sky = np.zeros_like(data).astype('float')
if verbose:
bar = ProgressBar(data.shape[0], ipython_widget=True)
for irx, row in enumerate(data):
fitted_sky[irx] = fit_chebyshev(row, degree, grow)
if verbose:
bar.update()
return fitted_sky
def fit_background(data, degree=5, grow=3, verbose=True, njobs=4):
"""
Parallelised background estimation for longslit CCD image
Params
------
data : array,
Input CCD data for background estimation.
degree : int,
Chebyshev Polynomial Order
grow : int,
Number of iterations to dilate around masked pixels
njobs : int
Number of processes to initiate for fitting
"""
kwargs={'degree': degree, 'grow': grow}
p = Pool(njobs)
fitted_sky = p.map(partial(fit_chebyshev, **kwargs), data)
return np.array(fitted_sky).astype('float')
In [3]:
ic1 = ImageFileCollection('../20170624/')
In [ ]:
blue_bias_list = []
for filename in ic1.files_filtered(obstype='Bias', isiarm='Blue arm', object='blue bias'):
print ic1.location + filename
ccd = CCDData.read(ic1.location + filename, unit = u.adu)
#ccd = ccdproc.create_deviation(ccd, gain=ccd.header['GAIN']*u.electron/u.adu,
# readnoise=ccd.header['READNOIS']*u.electron)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'])
blue_bias_list.append(ccd)
master_bias_blue = ccdproc.combine(blue_bias_list, method='median')
master_bias_blue.write('master_bias_blue.fits', clobber=True)
red_bias_list = []
for filename in ic1.files_filtered(obstype='Bias', isiarm='Red arm', object='red bias'):
print ic1.location + filename
ccd = CCDData.read(ic1.location + filename, unit = u.adu)
#ccd = ccdproc.create_deviation(ccd, gain=ccd.header['GAIN']*u.electron/u.adu,
# readnoise=ccd.header['READNOIS']*u.electron)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
red_bias_list.append(ccd)
master_bias_red = ccdproc.combine(red_bias_list, method='median')
master_bias_red.write('master_bias_red.fits', clobber=True)
In [ ]:
from astropy.convolution import convolve, Gaussian2DKernel
kernel = Gaussian2DKernel(25)
red_flat_list = []
for filename in ic1.files_filtered(obstype='Flat', isiarm='Red arm', object='well good flat r'):
ccd = CCDData.read(ic1.location + filename, unit = u.adu)
#ccd = ccdproc.create_deviation(ccd, gain=ccd.header['GAIN']*u.electron/u.adu,
# readnoise=ccd.header['READNOIS']*u.electron)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_red)
red_flat_list.append(ccd)
master_flat_red = ccdproc.combine(red_flat_list, method='median')
convolved_flat_red = convolve(master_flat_red.data, kernel, boundary='extend')
master_flat_red.write('master_flat_red.fits', clobber=True)
master_flat_red.data /= convolved_flat_red
master_flat_red.write('master_flat_norm_red.fits', clobber=True)
In [ ]:
blue_flat_list = []
for filename in ic1.files_filtered(obstype='Flat', isiarm='Blue arm', object='well good flat blu'):
ccd = CCDData.read(ic1.location + filename, unit = u.adu)
#ccd = ccdproc.create_deviation(ccd, gain=ccd.header['GAIN']*u.electron/u.adu,
# readnoise=ccd.header['READNOIS']*u.electron)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_blue)
blue_flat_list.append(ccd)
master_flat_blue = ccdproc.combine(blue_flat_list, method='median')
convolved_flat_blue = convolve(master_flat_blue.data, kernel, boundary='extend')
master_flat_blue.write('master_flat_blue.fits', clobber=True)
master_flat_blue.data /= convolved_flat_blue
master_flat_blue.write('master_flat_norm_blue.fits', clobber=True)
In [ ]:
ic1 = ImageFileCollection('../20170624/')
objects = ['USS422', 'USS337', 'USS7', 'SP1446+259']
for objname in objects:
print(objname)
"""
Blue Arm
"""
blue_target_list = []
for ifx, filename in enumerate(ic1.files_filtered(obstype='TARGET', isiarm='Blue arm', object=objname)):
print(ifx+1)
hdu = fits.open(ic1.location + filename)
ccd = CCDData(hdu[1].data, header=hdu[0].header+hdu[1].header, unit = u.adu)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_blue)
ccd = ccdproc.flat_correct(ccd, master_flat_blue)
ccd = ccdproc.cosmicray_lacosmic(ccd, sigclip=4.5, gain=ccd.header['GAIN'], readnoise=ccd.header['READNOIS'])
# Do sky subtraction
ccd.mask[:,785:800] = True
sky = fit_background(np.ma.array(ccd.data, mask=ccd.mask))
ccd.data -= sky
# Rotate Frame
ccd.data = ccd.data.T
ccd.mask = ccd.mask.T
blue_target_list.append(ccd)
#ccd.write('obj_'+filename, clobber=True)
blue_target = ccdproc.combine(blue_target_list, method='average')
blue_target.write('{0}_blue.fits'.format(blue_target_list[0].header['object']), clobber=True)
"""
Red Arm
"""
red_target_list = []
for ifx, filename in enumerateic1.files_filtered(obstype='TARGET', isiarm='Red arm', object=objname)):
print(ifx+1)
hdu = fits.open(ic1.location + filename)
ccd = CCDData(hdu[1].data, header=hdu[0].header+hdu[1].header, unit = u.adu)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_red)
ccd = ccdproc.flat_correct(ccd, master_flat_red)
ccd = ccdproc.cosmicray_lacosmic(ccd, sigclip=4.5, gain=ccd.header['GAIN'], readnoise=ccd.header['READNOIS'])
# Do sky subtraction
ccd.mask[:,785:800] = True
sky = fit_background(np.ma.array(ccd.data, mask=ccd.mask))
ccd.data -= sky
# Rotate Frame
ccd.data = ccd.data.T
ccd.mask = ccd.mask.T
#ccd.write('obj_'+filename, clobber=True)
red_target_list.append(ccd)
red_target = ccdproc.combine(red_target_list, method='average')
red_target.write('{0}_red.fits'.format(red_target_list[0].header['object']), clobber=True)
red_target.mask[785:800,:] = True
red_sky = fit_background(np.ma.array(red_target.data.T, mask=red_target.mask.T)).T
red_target.data -= red_sky
red_target.write('{0}_red_skysub.fits'.format(red_target_list[0].header['object']), clobber=True)
In [ ]:
for filename in ic1.files_filtered(obstype='Arc', isiarm='Blue arm', object='CuNe+CuAr b tar'):
hdu = fits.open(ic1.location + filename)
ccd = CCDData(hdu[1].data, header=hdu[0].header+hdu[1].header, unit = u.adu)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_blue)
ccd = ccdproc.flat_correct(ccd, master_flat_blue)
ccd = ccdproc.cosmicray_lacosmic(ccd, sigclip=4.5, gain=ccd.header['GAIN'], readnoise=ccd.header['READNOIS'])
ccd.data = ccd.data.T
ccd.mask = ccd.mask.T
ccd.write('arc_blue_'+filename, clobber=True)
for filename in ic1.files_filtered(obstype='Arc', isiarm='Red arm', object='CuNe+CuAr r tar'):
hdu = fits.open(ic1.location + filename)
ccd = CCDData(hdu[1].data, header=hdu[0].header+hdu[1].header, unit = u.adu)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_red)
ccd = ccdproc.flat_correct(ccd, master_flat_red)
ccd = ccdproc.cosmicray_lacosmic(ccd, sigclip=4.5, gain=ccd.header['GAIN'], readnoise=ccd.header['READNOIS'])
ccd.data = ccd.data.T
ccd.mask = ccd.mask.T
ccd.write('arc_red_'+filename, clobber=True)
In [ ]:
Fig, Ax = plt.subplots(2, 1, figsize=(12,8))
Ax[0].plot(np.arange(4099), np.median(CCDData.read('processed_data/arc_blue_r2544941.fit').data/convolved_flat_blue.T, axis=0))
Ax[0].set_xlim([1000,4000])
#Ax[0].set_ylim([0, 2000])
Ax[1].plot(np.arange(4096), np.median(CCDData.read('processed_data/arc_red_r2544942.fit').data/convolved_flat_red.T, axis=0))
Ax[1].set_xlim([1000,4000])
#Ax[1].set_ylim([0, 30000])
Fig.tight_layout()
blue_lines_pix = [1177, 1846, 1868, 1895,
2140, 2173, 2368, 2610,
2703, 2950, 3000, 3132,
3556, 3638]
blue_lines_w = [3273.96, 3858.58, 3868.53, 3891.98,
4103.91, 4131.72, 4300.1, 4510.73,
4589.9, 4806.02, 4847.81, 4965.08,
5330.78, 5400.56]
red_lines_pix = [1341, 1392, 1938, 1994,
2057, 2072, 2111, 2188,
2498, 2672, 2814, 3143,
3199, 3270, 3437, 3505]
red_lines_w = [5852.49, 5944.83, 6929.47, 7032.41,
7147.04, 7173.94, 7245.17, 7383.98,
7948.17, 8264.52, 8521.44, 9122.97,
9224.5, 9354.22, 9657.78, 9784.5]
In [ ]:
Fig, Ax = plt.subplots(1, 1, figsize=(6,6))
Ax.plot(red_lines_w, red_lines_pix, 'o')
Ax.plot(blue_lines_w, blue_lines_pix, 'o')